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Abstract 

Reduction of computational cost of solutions is a key issue to crack identification or crack propagation 
problems. One of the solution is to avoid re-meshing the domain when the crack position changes or when the 
crack extends. To avoid re-meshing, we propose a new finite element approach for the numerical simulation 
of discontinuities of displacements generated by cracks inside elastic media. The approach is based on a 
fictitious domain method originally developed for Dirichlet conditions for the Poisson problem and for the 
Stokes problem, which is adapted to the Neumann boundary conditions of crack problems. The crack is 
represented by level-set functions. Numerical tests are made with a mixed formulation to emphasize the 
accuracy of the method, as well as its robustness with respect to the geometry enforced by a stabilization 
technique. In particular an inf-sup condition is theoretically proven for the latter. A realistic simulation with 
a uniformly pressurized fracture inside a volcano is given for illustrating the applicability of the method. 

Keywords: Finite Element Methods, Fictitious domain methods, Xfem, Crack, Linear Elasticity Model. 
AMS subject classifications (2010): 74B05, 74S05, 65N30, 74R10. 


1 Introduction 

Recovering information on cracks located inside elastic media and predicting their propagation is a key topic 
for several geophysics and engineering problems. For instance, the analysis of ground deformation is used to 
search for the shape, position, and stress changes of magma-filled fractures at volcanoes (see [WCKdl2]) and 
of seismogenic faults [ONK+11] in order to assess the associated hazards. Studying the propagation of fluid 
filled fractures is central to the study of hydraulic fracturing of hydrocarbon reservoirs or the formation of ore 
deposits. 

These problems usually involve multiple computations for which only part of the boundaries is modified: De¬ 
termination of crack characteristics from ground deformation usually requires hundreds of computations with 
different crack configurations and study of fracture propagation requires computation of incrementally larger 
crack surfaces. 

For these problems to be numerically tractable, one possibility is to use methods in which the domain does not 
have to be meshed, such as boundary element methods [CMR83], for which meshing is limited to the boundaries 
so that modification of the system is only required for the modified boundaries. An other possibility is to use 
domain methods, such as fictitious domain methods in which meshes are non-conforming so that re-assembly 
of the whole system is avoided. 

Boundary element methods involve full non-symmetric matrices. These methods [CMR83] can take anisotropic 
media into account, but treating heterogeneous media requires the definition of new boundary, which increases 
matrices dimensions and involves very long computation times. As a consequence, heterogeneous media consid¬ 
ered using boundary element methods only have two different material properties. On the other hand, domain 
methods, such as finite element methods, can deal with anisotropic and heterogeneous media with no increase 
in the system dimension and no extra numerical cost. In order to avoid re-meshing when the crack is updated 
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and to treat heterogeneous media appropriately, we have chosen to develop a finite element method in which 
cracks are taken into account with a fictitious domain method. Besides accuracy, we require the method to be 
robust with respect to the problem geometry: The same results must be obtained whatever the way the crack 
is intersecting the mesh. 

Our method is based on an artificial extension of the considered crack (Figure 1) in order to split the domain 
into two sub-domains. It is different from the way the crack is extended in [BSES05], where the extension is made 
such that the crack withdraws into itself. The aim is to double properly the degrees of freedom around the crack. 
Discontinuities of the displacement field generated by traction forces (Neumann-type boundary conditions) im¬ 
posed on both sides of the crack have to be simulated. We can impose other Neumann-type conditions on this 
interface. The fictitious domain approach we implement is inspired by Xfem [MDB99], since it consists partially 
in cutting the basis functions near or around the interfaces; But, unlike Xfem [CLR08, LPRS05, LYMS14, FB10], 
we do not enrich our finite element spaces functions with singular functions, as we intend to avoid these en¬ 
richments by minimizing the number of updates when the position of the crack is modified. This approach 
has been first introduced in [HR09] for the Poisson problem, and in [CFL14] in the context of Fluid-Structure 
Interactions, for Dirichlet conditions. 

Across the artificial extension considered for the crack, there is no discontinuity, so we impose a homoge¬ 
neous displacement jump condition with a Lagrange multiplier. This boundary is also taken into account with 
a fictitious domain approach, and the multiplier aforementioned is a dual variable for which we want a good 
approximation. This point is crucial for two reasons: First this quantity represents the constraints at stake in 
this region, whose the knowledge is required for crack propagation problems. Secondly this quantity takes part 
in some inversion algorithms (like gradient algorithms based on the computation of a direct adjoint system). 
Getting a good approximation of this multiplier is not guaranteed a priori , because the degrees of freedom 
considered on the underlying mesh do not match the crack and its extension, and it is the price to pay when 
we do not provide enrichment of basis functions like it is done with Xfem-type methods, where the enrichment 
of the basis functions has to be done according to the geometry of the crack, which makes part of the up¬ 
dates we want to avoid. For circumventing this drawback, we carry out a stabilization technique of augmented 
Lagrangian-type a la Barbosa-Hughes [BH91, BH92]. This technique theoretically ensures an unconditional 
optimal convergence for the multiplier, but improvement is mainly observed when we perform numerical tests 
related to the robustness with respect to the geometry. 

The paper is organized as follows: In Section 2 we introduce the theoretical problem, we explain why and how 
we adopt the extension method for considering the crack and we give the continuous abstract weak formulation. 
In Section 3 we first detail the discrete formulation, in particular the principles of the fictitious domain method 
developed for the interfaces associated with the crack. Then, in Section 3.2, we provide the theoretical analysis 
without the augmented Lagrangian technique. Next in Section 3.3 we introduce the stabilization technique which 
forces the multiplier to reach the desired value. Lemma 3 shows that an Inf-Sup condition is automatically 
satisfied while performing the stabilization technique. Section 4 is devoted to practical explanations for the 
implementation. In Section 5.1 we present some 2D numerical tests without stabilization which estimate the 
rates of convergence for the displacement as well as for the multiplier. First illustrations in 2D are also given. 
The choice of the stabilization parameter is discussed in Section 5.2, and this part is concluded with numerical 
tests in Section 5.3 providing rates of convergence with the stabilization. In Section 6 we show that the interest 
of the stabilization technique lies in the criteria of robustness with respect to the geometry. Last, a realistic 
3D simulation of an over-pressured magma-filled fracture inside Piton de la Fournaise volcano is performed in 
Section 7. Conclusion is given in Section 8. 
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2 Setting of the problem 

2.1 The original elastic problem 

Given a domain of (d = 2 or 3), and a crack IV CC 0 represented by a curve or a surface parameterized 
by an injective mapping, we consider a static linear elasticity model governed by the following system: 

{ — div <jl{ u) = f in Q, 

u = 0 on dQ, 

01 ,( 11)11 =pn on Ty. 

In this system the displacement of the solid is denoted by u, body forces (like the gravity) by f, and 0z,(u) = 
2hls(vl) + A/,(div u)I R d denotes the Lame stress tensor, with 

e(u) = i(Vu + Vu T ). 

The real numbers /il and A l are the Lame coefficients. The pressure force of value p > 0 is applied on both 
sides of the crack Tt, so we have to specify the outward normal n on Tt- In order to determine solutions on 
both sides of the crack, we relate displacements on each side of the fracture to the displacement discontinuity 
across the fracture. 

2.2 Extension of the fracture 

In order to give a sense to both sides of the crack, we have to be able to determine whether a point of the domain 
lies on one side of the fracture or the other. The most convenient way we have found consists in uncoupling 
the problem by setting two unknowns displacements instead of a global one. We extend the crack Ft to T, as 
shown in Figure 1 below: 



Figure 1: Splitting of the domain according to the crack. 

Let us justify why we adopt this method of extension of the crack. First, in view of the boundary conditions 
we impose on the crack Tt, it is necessary to define two unit normal vectors on this interface, and thus to be able 
to split the domain as a function of the discontinuity induced by these boundary conditions. As aforementioned 
in the introduction, it is then more convenient to consider a configuration for which we are able to locate a 
point with respect to the crack. 

The relevance of such an approach can be underlined by the theoretical context. Indeed, the Cauchy problem 
as well as the identification for crack-like boundaries is originally ill-posed. See for instance [StaOl]. Splitting 
the domain as described above can enable us to recover existence and uniqueness of a (at least weak) solution. 
The question of the possibility of making a crack extension has to be answered. Namely we seek the condi¬ 
tions under which a crack extension splitting domain Fl into two sub-domains exists. A first answer could be 
to only consider cracks for which an extension exists. Such cracks necessarily have to be represented as an 
injective curve/surface. An argument for this answer is more physical and results from the fact that, modeling 
crustal displacements, the domain considered includes cracks. In [CDF05] fractures inside a realistic 3D volcano 
were described by a specific parameterization, which can easily be extended. Moreover, this extension can be 
constructed and implemented in practice. For more details on this kind of extension, we refer to Section 7.1. 
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2.3 The transformed problem 

The global domain Q is now split into two sub-domains and . We have: 

r = r 0 n r T , n = n + uru 


Let us now denote u + = U|^+ and u = U|^+. On the artificial boundary To - which is not connected - we 
have to ensure the continuity of the displacement, namely u + — u - =0. The system becomes: 1 


—div < jl ( u ± ) = f 
u ± = 0 
(cr L (u)n) ± = pn ± 

M=0 

Mu)]n+ = 0 


in 

on dfi n d 
on Tt, 

across To = T \ Tp, 
across T 0 . 


(1) 


Notation [ip\ = <p + — refers to the jump of a function cp across Tq. We assume that f G L 2 (D), and we 
still denote by f its restriction to or Yl~. The homogeneous Dirichlet condition on dYl can be replaced by 
non-homogeneous boundary conditions mixing Neumann conditions and Dirichlet conditions. Moreover, the 
pressure condition on Yt can be replaced by general conditions, as (cri / (u)n) ± = ±g. 


2.4 Continuous formulation 

Consider the following functional spaces: 

V+ = {v G H 1 (D + ) | v = 0 on <9D H <9D+} , 

V- = {vGtf(O-) | v = 0 on <9Dn<9D - } , 

w = (h 1 / 2 ^))'. 

We choose to impose the jump condition on Tq by a Lagrange multiplier A. A weak solution of system (1) can 
be seen as the stationary point in V + x V - x W of the following Lagrangian: 

£ 0 (u + ,u“,A) = \( <jl(u + ) : e(u + )dfi + + i / cr L (u“) : e(u _ )dfT 

1 Jn+ 1 Jn- 

— f f • u + dD + — / f • u _ dD _ — f u + pn + dT T — f u _ -pn _ dT T 

Jn+ Jn- jvt Jrv 


+ (A,(u+ - u )) 


W,W' 


( 2 ) 


In this expression ( • , • )w,W' denotes the duality pairing between W and W', and : s(u) = 

trace (<ji / (u)5(u) t ) denotes the classical inner product for matrices. Let us recall that the bilinear form 
(u, v) i-G u) : e:(v) is symmetric. 

Remark 1. Note that in the writing of this Lagrangian the jump condition [<tl(u)] n + =0 onY o (fifth equation 
of (1 )) is no longer taken into account. Indeed, the first-order optimality conditions for Co give 


f <7z,( u± ) : e(v)dfi ± ±(A,v) ww , = f f • vdf2 ± + f ■ 

Jn± ’ J n± Jv T 


v•pn^dT^, 


for all test function v G V^. On the other hand, taking the inner product by v of the first equation of (1) 
yields, after integration by parts 


j 

Jn- 


Thus we identify A 


) : e(v)dD ± - (a L ( u ± )n ± , v) w?w/ 
-(Jl{ u + )n + and X = aL(u~)n~. 


f • vdfT 




J r T 


v • prCdr^. 


x The symbol =b represents the fact that we consider both formulations involving the symbols + and —, for the sake of concision. 
The outward normal of domain is denoted by nA 
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The variational problem derived from the functional as first-order optimality conditions is then 


where we set 



Find (u + , u , 

A) in 

V+ x V- 

x W such that 


f ^((u + ,u 

~i A); 

v) = Z±(v) 

Vv e V ± , 


1 Bo((u + ,u- 

, A );n) = 0, 

V/i € w. 

*4o ((u+ 

c 

1 

V 

II 

/ 

ctl( u *): £ 

(v)dfi ± ± (A,v) W: 



Jn± 



^(v) = 

/ 

f •vdO ± + 

/ v•pn ± drT, 


Jn± 


Jr T 

#o((u + , 

u“,A );/x) = 


u+\ — 

/W,W' 



(3) 


(4) 


3 Discrete formulation 

The discrete formulation is adapted from [HR09] and [CFL14]: It is a fictitious domain method, in which the 
choice of the degrees of freedom for the multiplier on the boundary To is made independently of the mesh. Let 
us first explain how we proceed to take into account degrees of freedom which do not originally he on the edges 
of the mesh. 


3.1 The fictitious domain approach 

Unknowns for the fictitious domains are first considered on the whole domain Q. Let us consider some discrete 
finite element spaces, C H 1 (U) and C L 2 (U). These spaces can be defined on the same structured 
mesh of U, which can be chosen Cartesian. We set 

Vh = {v h ec(fi) | v h]&a = 0, v h{T e P(T), VTe%}, (5) 

where P(T) is a finite dimensional space of regular functions, including polynomial functions of order k > 1 on 

a triangle T in the set Th of the triangles of the mesh. See [EG04] for more details. The mesh parameter stands 

for h = max hr , where hr is the diameter of the triangle T. We define 
TeT h 


V£:=V h |n + , V,, :=V/,|o-, W h := W h]Fo , 

which are natural respectively discretizations of V + , V - and (H 1 / 2 (T 0 )) / . The selection of degrees of freedom 
for these spaces is illustrated in Figure 2. This approach is similar to the extended Finite Element Method 
[MDB99], except that the standard basis functions near the boundary T are not enriched by singular functions 
but only multiplied by Heaviside functions (iL(x) = 1 for x G and iL(x) = 0 for x G 9 \ and the 
corresponding products come up in the integrals of the variational formulation of the problem. This kind of 
strategy is also adopted in [GW08] and [CHM12] for instance. 
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Figure 2: Selection of degrees of freedom: The black base nodes are kept for displacements u^, the red ones are 
used for cutting the standard basis functions, and the yellow ones are those of the multiplier A. 


An approximation of problem (3) is given as follows 

Find (u£, u^, A h) in x x such that 

Vv+ G V+ 


«h u b v P+ b o( x k, vp=z+(vp 

«o K> V J -&o(Ah,v h ) =l-(v h ) 
b o(u£,v h ) - b 0 (u h ,n h ) = 0, 


Vv^GV^, 

V/*/, € W h , 


where 


J o( u+ > v+ ) = / a L (u+) : e(v + )dft+ a 0 (u ,v )= [ a L (u ) : e(v )dfi , 

Jn+ Jn- 


bo(u ± ,n h )= , 
Jr 0 


/x • u ± dr 0 . 


(6) 

(7) 

(8) 


Note that the duality pairing between (H 1 / 2 (Fo)) / and H 1 / 2 (Fo) has been turned in the expression of bo into the 
inner product of L 2 (T 0 ). We could avoid this by using a Laplace-Beltrami operator, but under strong regularity 
assumptions, it is simpler to proceed like this. Thus we now consider that 

w' = w = L 2 (r 0 ), w h c L 2 (r 0 ), 

but we can keep the abstract formalism between W and W' in the mathematical analysis. 


3.2 Theoretical convergences - limitations on the orders 

We make the following hypothesis: 

(H): JI h €W h : 6o(v ± ,Mj=0Vv ± eV± =* Ji h = 0. 

This hypothesis is not as strong as an Inf-sup condition which would lead automatically to the optimal order 
of convergence for all the variables - primal and dual. See [BF91] for more details. It only requires that the 
discretization space for the displacement is at least as rich as the one chosen for the multiplier on To. 

Lemma 1. The bilinear forms a,Q and defined in (7) as 

a± : (u, v) i y / (u) : 5(v)dQ ± 

Jn± 

are respectively uniformly -elliptic and -elliptic. Namely there exists a > 0 independent of h such that 
for all v± G we have 

a o( v± > v± ) > a || v dlv± • 
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Proof. Since C V ± , let us show the coerciveness of a ± in V ± . Arguing by contradiction, assume that there 
exists (v±) such that for all n G N we have 

n[ ol(v±) : e(v^)dn ± < ||v,‘, ||^ ± . (9) 

J n± 

After a quick calculation we have 

o-l(v^) : e(v±) = 2^ L |e(v±)| 2 <iXti + A L (div v±) 2 . 

On the other hand, because of the homogeneity of the inequality (9), we can assume that ||v±|| v± = 1, without 
loss of generality. Then, inequality (9) implies in particular that £(v±) tends to 0 in [L 2 (^ ± )] dxd . The Rellich’s 
theorem gives us a subsequence (v±) which has a limit in L 2 (f] ± ). From the Korn’s second inequality (see 
for instance [EG04], Theorem 3.78, or [Cia97, page 10] for the proof) there exists C > 0 such that 2 

R-^|| Hl(n±) - C (ll v m-vt|| L2(n±) + |k(v^)-e(v±)|| [L2(n±)]d><d ). 

Thus (v^) m is a Cauchy sequence in H 1 (fl ± ), and so it tends to some which satisfies in particular s(v±) = 0. 
It is then well-known that reduces to a rigid displacement, that is to say is an affine function of type 
(x) = ^ + r ± A x (see [Tem83, page 18]). Besides, the trace theorem implies that we have also = 0 
on It enables us to conclude that ^ = r ± = 0 and thus = 0 in which belies the assumption 

ll V n|lv± =1 - D 

Proposition 1. Assume that the hypothesis (H) holds. Then Problem (6) admits a unique solution (u^, u^, Xyf). 

Proof. Problem (6) is of finite dimension, so the existence of the solution is a consequence of its uniqueness. For 
proving uniqueness, it is sufficient to show that (u£, u^, A/J = 0 when Z ± = 0. If we add the first two equations 
of (6) while choosing and = u^, the third equation of (6) enables us to write 

«o ( u £> u £) + a o ( u h> u /0 = 0, 

and thus by Lemma 1 we get = 0. Then, the first two equations of (6) reduce to the conditions of 

Hypothesis (H), and thus we can conclude that \ h = 0. □ 

We now define the space 

< = {(V+ v k -)eV+xV k - |6o(vt-v^ )/th )=0V/* h 6W h }. 

Let us give the abstract error estimate for the displacement. 

Proposition 2. Assume that the hypothesis (H) is satisfied. Let (u + ,u _ , A) and (u^,u^, A/J be the respective 
solutions of Problems (3) and (6). There exists a constant C > 0 - independent of h - such that 


u — u 


h |lv+ 


+ u — u 


h IIV- 


< 


C inf 


u+ - V /Mlv+ + 


u — V 


dlv-) + J|wJ |A ^ llw ) 


( 10 ) 


Proof. From Lemma 1, for all test functions (v^, ) G we have 

a ll u ft _ v dlv± ^ - v b u h - v h) 

= A(. u± -Vh,U-h - + K U h - ^h) - a ?( u± ,Uh - v^) 

= 4^ - A A - A) ± - v±> w w ,. 

Summing the two inequalities above which correspond to the symbols + and —, and denoting — v^, 

we get 

a (ll w dli+ + ll w *lli-) ^ a+(ut^v+ w+) + ao(u--Vfc,w^) 

+ <A,(u+-u^)-(v+-v^)) Wi w'- ( n ) 

2 In the following, C denotes a generic positive constant independent of the mesh size h. 
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From the definition of , and noticing that (u£, ) eV°, we can write 


(A.fu,;-- u,,) ( V/ ; v /( )) w w; = (\-n h ,(u+ -U h ) - (v+ - vj) 


W,W' 


Thus, for (v+, ) e Y° h and n h e W h , the inequality (11) becomes 


a w 


h ||V+ 

2 


+ W' 


w 


h Ilv+ 


+ W 


h 11V — 
2 


) < 4(u+ - v+,w+) + a 0 (u -v h ,w h ) + (A-Ai A , 


w/ 


v» h e w h . 


W h )w,W' 


h IIV- — 


< C ( u 4 


h Ilv+ 


U 


h IIV — 


+ ||A - AiJ w ) (||w+|| v+ + ||w h || v _) , 


and, recalling that w^ 


Uu 


v±, we obtain (10). 


□ 


These results show us that, under Hypothesis (H), Problem (6) admits a unique solution, which satisfies the 
a priori estimate (10). We have no such estimate for the approximated multiplier A^. 

Besides, the estimation of the convergence rate for the primal unknown given in [HR09] for the Poisson problem 
can be transposed to our elasticity problem. Proposition 3 of [HR09] - page 1480 - gives an order of convergence 
at least equal to 1/2. It can be adapted to our case as follows. 

Proposition 3. Assume Hypothesis (H). Let (u + ,u _ ,A) be the solution of Problem (3) such that u ± G 
H^/2+i+£(q±) p| y± j or some £ > o. Assume that 


inf ||A-mJw ^ Ch& 

h 


for some S > 1/2. Then we have 


ll u±_u ^llv± - C ^ 

Proof As it is shown in the proof of Proposition 3 of [HR09], page 1481, for any u ± G H d / 2+1+£ (Q ± ) D V 1 ^, 
we can show that there exists two finite element interpolating functions (v£, v^) G such that 


u± - v h\\v± ^ c ^- 


( 12 ) 


Such functions can be constructed as standard interpolating vectors of (1 — 7^/ i )u ± , where 77 p is a cut-off 
function equal to 1 in a vicinity of the boundary To. More precisely, in a band of width 3/i/2, vanishes on 

all the convexes intersected by To. So we can be sure that vanishes on Tq, and that (v£, v^) G . The 

announced estimate can then be deduced from (10) combined with (12) and the hypothesis on the interpolating 
functions fi h in the statement of this proposition. □ 

This theoretical limitation of the order of convergence is very common for fictitious domain methods (without 
special treatments as those provided by stabilization techniques). However, the numerical tests provided in 
Section 5.1 indicate that this rate of convergence is not sharp enough; Indeed, a test which would show the 
optimality of this estimate is hard to provide: It has been provided in [HR09] in a very specific configuration, 
and merely not found in [CFL14]. In our case, we are not able to observe numerically this limitation when the 
outer boundary dQ fits the mesh, but we observe it easily when this boundary and the Dirichlet condition we 
impose on it are taken into account with a fictitious domain method in addition to the conditions considered 
on T. We do not comment furthermore this technical point, since it is beyond the scope of this paper. 


3.3 Stabilization technique 

The augmented Lagrangian technique 

In order to enforce convergence on the multiplier A for all fracture configurations, we develop a stabilization 
method a la Barbosa-Hughes [BH91, BH92]; The idea is to use the Lagrangian functional Co given in (2) with 
an additional term which takes into account the required value for A. Thus we now consider the Lagrangian 
functional below 

£(u + ,u - ,A) = £ 0 (u + ,u“, A) - j- [ |A + cr L (u + )n + | 2 dr 0 — ^ [ \A - cr L (u - )n - | 2 dT 0 . (13) 

z ^r 0 z J r 0 
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Observe that this extended Lagrangian is equal to the previous one for an exact solution. The additional 
quadratic term tends to force A to reach the wanted value of the surface force ox( u - )n - = — ox( u + )n + (see 
Remark 1). The parameter 7 > 0 indicates the importance given to this requirement. This additional term 
deteriorates the positivity of the Lagrangian £, and thus the coerciveness of the formulation which stems from 
it. As we cannot choose a too large 7 , this technique is not a penalization method. This choice is discussed in 
Section 5.2. 


A formal calculation of the first variations gives 




(v) = -7— r(v) — 7 / A • ai(v)n + dr 0 — 7 [ o- L (u + )n + • o- L (v)n + dr, 

^ u+ Jto J r 0 

( v ) + 7 / A • (JL(v)n _ dr 0 - 7 / (JL(u _ )n _ • (JL(v)n _ dr l 
JTn J r n 


II 

<5£o 

du+ 

sc , , 

SCq 

II 

1 

3 

<5u _ 

sc, , 

SC 0 

m'" 1 = 

SX 


^(v) -27 A • MdlV 

The stabilized discrete formulation is then: 

Find (u^,u^,Aft) in x x Wft such that 

f •^ ± (( u t> u h>A/ l );v /l ) = l ± (v h ) Vv h e V±, 

\ =0, \//j, h eW h , 


where 


. 4 . 1 ({n ‘. u , A): v) 


A•vdT 0 


[ cr L (u ± ) : 5(v)df] ± ± [ 

J ry 

T7 / A • cr L (v)n ± dr 0 - 7 / cr L (u ± )n ± • cr L (v)n 
^r 0 J r 0 

#((u + ,u _ , A); /x) = [ ix u + dr 0 -/ /x - u _ dT 0 

Jrn 7r 0 


o, 


'r 0 

-7 


1 */10 

/ /i • a L (u + )n+dr 0 + 7 f /it • ( 7 L (u _ )n _ dr 0 - 27 /* AjudlV 

Aq 


(14) 


The discrete Inf-Sup condition 

Let us choose 7 = 70 h, where 70 > 0 is independent of h. Note that Problem (14) can be rewritten as follows: 
Find (u^,u^, Aft) G x x Wft such that 

■W(K,u^A k );(v+,v ft -,^)) = V(v+,v^,^) € V+ xV( x W h , 

where 

Al((u + ,u _ , A); (v + , v _ ,/i)) = f cr L (u + ) : 5(v + )dfl + + [ a L (u~) : s(-v~)dQ~ 

Jn+ Jn- 

+ [ A (v + - v _ )dr 0 + [ ix - (u + - u“)dr 0 
Jr 0 J r 0 

- 7 oh / (A + cr L (u + )n + ) • (/x + cr L (v + )n + )dr 0 - 7 0 h / (A - cr L (u _ )n _ ) • (/x - <T L (v _ )n _ )dr 0 , 

Jr 0 Jr 0 

7 r (v + , v _ , /x) = f f • v + dfl + + f f • v _ dfl _ + f v + -pn + drT+ [ v~-pn~dFT- 

Jn+ Jn- J r T J r T 

In the following, we will need an assumption for the theoretical result stated in Lemma 3: 

(A): /i||cr L (v±)n ± ||L 2(ro) < C||v±||^ ± Vv± e V±. 

This assumption is of same type as the ones considered in [HR09] (page 1482, the equation (27)) for the Poisson 
problem, or in [CFL14] (page 84, Assumption A1 stated in Section 4.2) for the Stokes problem, both for the 
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theoretical study of the fictitious domain approach stabilized a la Barbosa-Hughes. 


We will also use an estimate for the L 2 -orthogonal projection from H 1 / 2 (Fo) on W^. It has been established 
and proved in [CFL14] (page 85, Lemma 3): 

Lemma 2 . For all v G H 1 / 2 (ro) we have 


ll^ v - v ll L 2(r 0 ) < CVh ||v|| Hl/ 2 (ro) , 
where Vh denotes the L 2 -orthogonal projection from H 1 / 2 (Fo) on W^. 

The main purpose of this subsection is to prove the following inf-sup condition. 

Lemma 3. Assume that (A) holds. For 70 small enough, there exists C > 0 independent of h such that: 


inf 


sup 


( u ^ u fc> A 0 GV ii xV /* xW h (v+v-,^)ev+xV-xw h Ill u ^ u /i j^IH \\\ v hi v hi Vhl 
where the norm 111 • 111 is defined as follows: 




> c. 


u + ,u ,A| 


= u 1 


2 

v+ 


— 112 


1 


„v-nH U U ~ 


— 112 


lL 2 (r 0 ) 


L 2 (r 0 ) 


Proof Ideas for this proof are similar to the ones used in [HR09] (page 1482, Lemma 3) or [CFL14] (page 85, 
Lemma 4), and inspired by [Ste95] (page 144, Lemma 6). First, we estimate 


M (( U h’ U h’ X h);(ul,u h ,-X h )) = Go ( U £> U fe) + «0 ( U h> U J + 270 h 

J rn 


'dr 0 


-70 h \a L (u^)n + \ 2 dTo-joh \a L (u h )n 

Jr 0 J r 0 

(Il u h|lv+ + ll u h||v-) +2 7o/i f |A h | 2 dr 0 

JFq 


‘dT 0 


> a u 


-C'l'o (||u+||L + ||uj|v-) 

> f (|| u h|lv+ + || u fc|lv-) + 2 ^ h H A ^llL 2 ( r 0 ) - 

by using Lemma 1 and Assumption (A), and then by choosing 70 small enough. Next, let us consider ~p h = 
pV h (u+ — u^). In view of Assumption (A), the Cauchy-Schwarz and Young inequalities, we have 

(0,0,/Z h )) = t \\p h (-u+ -u-)||^ 2(ro) -7o/i^ cr L (u + )n + • JZ h dT 0 


+ 7 o h cf l {u )n • M/idr 0 - 270ft / X h -^ h dT 0 

J r 0 J r 0 

- \ \\ Vh ( u h ~ u ^)llL2(r 0 ) “ 270 ^ IIV|| L 2(r 0 ) ||/*/,|| L 2(ro) 

-7o (\/ft||az,(u+)n + || L2(ro) +\Zft||cr L (u-)n-|| L2(ro) ) \/ft||^|| L2(ro) 

^ \ \\ V h( u h - u ft)|lL (ro ) -7oft||A ft ||L2 (ro) 

II V hK - u^)||^ 2(ro) - C70 (||u+||^ + + ||Ufc ||v_) 

> ^ \\ Vh ( u + - u ^)\\ l 2(ro) - J 0 h \\\ h \\ l 2(ro) - C^o (||u+|| 7 + + Ibhllv-) * 

for 70 small enough. Now, choosing (v^, v^, /j, h ) = (u^, u^, — A^ + V~Fh) ~ where 77 > 0 is supposed to be small 
enough - the previous inequalities above yield 

M( u h>“fc> A /0;( v h> v fc >**/,)) > ^ ||^( u 4 — u ^)|lL2(r 0 ) + ( 2 7oft-??7oft) ||A h ||L 2 ( r 0 ) 

+ (| -Crr/o) (||ufc||v+ + IKIIv-) • ( 15 ) 
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The projection term can be handled by using Lemma (2) as follows: 


\Vh{n 


) 11 L 2 (F 0 ) — 


I v h{n+-n h ) 


h blL 2 (To) 


= u; -u 


> 

> 


- u 


u: 


u 


h 11L 2 (F 0 ) 
2 

h 11L 2 (F 0 ) 
2 

h 11L 2 (F 0 ) 


- \\ v h(u+ -uj-(u+ -u 


-■) 2 

h> L 2 (r 0 ) 


-C/i||u+-u^|| H 1 / 2 (ro) 


■ 2 Ch ( u 


x h Ilv+ 


U 


h 11V — 


(16) 


Then, in view of (15) and (16), we have 

> (| - Crno - 2?yc) ([|u; 


+? 


h Ilv+ ll U h llv-) 

l IK - u h IIL 2 ( r o) + K 2 - v)h ||A fc ||L 2 (ro) , 

and thus, recalling that (v£, v^, ii h ) = (u^, u^, —A h + V~Ph)i an d choosing 70 > 0 and 7 > 0 small enough we 
obtain 


M((u +, u h , X h ); (v£ , , n h )) > C 11 |u+, U/l , A, 

where C > 0 is independent of h. On the other hand, we have 




< M\\\u+,u h ,X h 


(17) 

(18) 


where M > 0 is independent of h too. Indeed, 


v h v h»/*/, 


< 


< 


l u ^ u /^ A /*lll + 7 7lll 0 >°,M/JII 

u / ;.u /( .A,,||| + \\Ph{n+ - Ufc)|| L2(ro) 


K U /,.>V||| + -K IK -u 


»? 


< (l+»7) K U K 


\[h 


^ IIL 2 (r 0 ) 


Thus, combining (17) and (18) gives the inequality 






C 

M 


> —UK, u ft ,AJ 


which leads to the desired result. 


□ 


Note that the bilinear form A4 is bounded for the triple norm on V + x V - x W uniformly with respect to 
the mesh size h. Then Lemma 3 leads to the following error estimate with a Cea type lemma (see [EG04] for 
instance): 


u + - U+U — U, , A — AJ 


< C inf 

(^v-,/x,j€V+xV-xWft 


u+ -+> u - v h. A -^I 


From that we can invoke the extension theorem for the Sobolev spaces, the standard estimates for the nodal 
finite element interpolation operators and the trace inequality 

1 , 


II v ||l 2 (Fo) < C (^IM| L 2 (T) + ^||v||l 2 (T) 

on any convex T G Th (see Appendix A of [HR09], page 1496, for more details). As a consequence, it yields the 
following abstract error estimate 


max( || u — u 

< 


c(h k “ ||u 


h Ilv+ * 
+ 11 


u — u 


+ 1 


lH fc n + l(0+) 


h Hv+ 


\ _ \ + 

5 h L 2 (r 0 ) 


lH fc u+i(0-) 


) < |||u 4 

■ h kx 1 II. 


u h ,X-X h 


llH fc A+i(r 0 ) 


)• 


where k u and k\ are the respective degrees of standard finite elements used for the displacements u^, and 
the multiplier A. 
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4 Practical remarks on the implementation 

4.1 Matrix formulations 

In matrix notation, the formulation (6) gives 


( A + 

0 

Bf 

0 

1 o 

~Bf 

\w~ 


0 



F+ 

F" 

0 


where U + , U and A are the degrees of freedom of u£, and A h respectively. The matrices Aq , A 0 , Bq 
and Bq are the discretization of (7)-(8). If we denote the basis functions of the spaces V^, and by 
{< p \}, {<£> i } and {V’i} respectively, we have: 


= j Q+ vdvt) ■ e(v>+)dfl+ 

Wb = / r 


i • <Pj dr o, 


(4) )ij = J n _ a L(Vi ) : £ (Vj ) dfr 


i dr 0 . 


Vectors F ± are the discretization of (4). 

In matrix writing, the stabilized formulation (14) corresponds to 


' 4+ 

0 

b +t 

0 

A~ 

-b~ t 

B+ 

-B+ 

-c 



F+ 

F“ 

~o~ 


(19) 


The matrices A , B± and C introduced here are the respective discretizations of the following bilinear forms: 

a ± (u ± ,v) = [ cr L (u ± ) : 6 (v)df] ± - 7 / <r L (u ± )n ± • cr L (v)n ± dr 0 , ( 20 ) 

Jn± J r 0 


6(v, A) 


[ A • vdT 0 =F 7 [ 
-T 0 J r 0 


A • vdT 0 =F 7 / ^ cr L (v)n ± dr 0 , 


c(A,/x) = 2 7 [ . 

-To 


A • /xdTo. 


( 21 ) 

( 22 ) 


4.2 Solving the whole system 

Let us recall the extended stabilized formulation introduced in Section 3.3: 

Find (u + ,u“, A) in V + x V - x W such that 
( a+(u+, v+) = — 6 (A,v+) + Z+( v+ ) Vv+ e V+, 

< a“(u“, v“) = b(X, v _ ) + Vv" g V“, (23) 

[ b( u+, fi) - b(u~, fi) - c( A, ii) =0, G W, 

where the linear forms are given by (4) and the bilinear forms a^, b and c are given by (20)-(22). Recall also 
that in the discrete formulation (6) we considered for the bilinear form (8) the inner product in L 2 (Tq) instead 
of the duality pairing ( • , • )w,w / ? leading to consider that W = W' = L 2 (To). 

In view of the formulation (23), we deduce that the mapping A i— u ± (A) is affine, and u ± (A +tji) = u ± (A)+tu; ± , 
where G are such that 


a + (u> + , v + ) = — 6(/Z, v + ) Vv + G V + , 

, v _ ) = fe(/Z, v _ ) Vv - G V - , 

c(/Z, /i) = 0 V/iG W. 

If we set v 1 ^ = u ± in the two first equations of (23), we get after summation 

a + (u + , u + ) + a - (u - , u _ ) = — 6(A, u + ) + 6(A, u _ ) + / + (u + ) + / - (u - ). 


(24) 

(25) 
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Substituting this equality in the expressions of the Lagrangian (13) based on the expression (2), we obtain the 
following dual functional 

J*(A) = £(u+(A),u-(A),A) 

= — ^ a+ ( u+ > u+ ) - u - ) - ^c(A, A). 

The dual problem is therefore the maximization problem below: 

Find A* G W such that: J*(A*) > J*(/x), V/xGW. (26) 

For a given direction ju, the directional derivative for J* is given by 

DJ*(X).Ji = — a + (u + , tj + ) — a _ (u _ , u> _ ). 

In view of the sensitivity system (24)-(25) with = u^, and in view of Remark 1, this derivative can be 
expressed as 

= b + (u + ,Jl) — b~(u~,JI) = (ji, (u + — u“)) w w , . 

Thus we deduce that the gradient of J* - for the L 2 (Fo)-inner product - is DJ*(\) = [u]. Moreover, if JI is a 
search direction for J*, we obtain 

DJ (A = (/x, [u] + t [u?]) w ^w'• 

The optimal step-size t* in the search direction JI is therefore 

£* (Ab [ u ])w,W / 

(/X, [<*>]) w?w , 

Since J* is quadratic and (strongly) concave, a good algorithm for solving (26) is a conjugate gradient algorithm. 
The Fletcher-Reeves conjugate gradient algorithm for the maximization of the dual functional (26) is summarized 
in Algorithm 1. 


Algorithm 1 Uzawa conjugate gradient algorithm 


Initialization: k = 0. For X 0 given, compute uo such that 

« + ( u o - v+ ) = -b(v + ,\ 0 ) + l + (v + ) 

Vv+ G V+, 

a~(uo,v~) = A (l ) — / (v ) 

Vv" € V". 

Set g 0 = [uo] on To (gradient), and JI 0 = g 0 (search direction). 

Iteration k > 0. Assume that u kl g fc , Ji k are known. 

Sensitivity: Compute such that 

a+ ( u, fe> v+ ) = -H v+ ,74) 

Vv+ G V+, 

a ”( w fe> v ”) = & ( V “-Pfc) 

Vv“ G V“ 

Step size: t k = [u fe ]) L 2 (ro) /(p fc , [w fc ]) L 2 (ro) 

Update: u± +1 = u± + t k w£ 

Gradient: g fe+1 = g fc + t k [u k ] 

Direction: JI k+1 = g k+1 + /3 k JI k , where /3 k = (g fe+1 , g fe+1 ) L 2(r 0 )/(g fc , gfck^ro) 


We iterate until the norm of the gradient becomes sufficiently small, namely 

(S/c5 S/c)l 2 (Fo) ^ ^ §o)l 2 (To) ’ 

for some tolerance 5 > 0. The main characteristic of the dual algorithm is that, at each iteration, we solve one 
uncoupled problem. 
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4.3 Libraries used in the implementation 

Numerical implementation is done with the code developed under the Getfem++ Library (see [RP]). As 
mentioned previously, our program is based on the approach initially introduced for the Poisson problem (see 
[HR09]). In dimension 2, solving the global system can be made by using the library SuperLU [DGL]; In 
dimension 3 the Uzawa iterative algorithm given in Section 4.2 and adapted for system of type (19) is carried 
out while using the GmmH— b Library (included into the GetfemH— b library). 

For this kind of boundary value problems, in Getfem++ several difficulties have been tackled. These include: 

- Defining basis functions of from traces on To of the basis functions of W^. In particular, the 
independence of these functions is not automatically ensured, and redundant degrees of freedom have to 
be eliminated to avoid handling non-invertible systems. 

- Localizing the interfaces corresponding to the crack and its extension. For that purpose, a level-set 
function method is implemented in the library. 

- Accurately computing the integrals over elements at these interfaces during the assembly procedure, this 
involves calling the Qhull Library (see [BDH96]). 

4.4 Tools for defining and updating the geometric configuration 

In practice, cracks are numerically described with level-set functions. A first one Is i describes the whole bound¬ 
ary T = Ft U To with the equation Isi = 0. Then, we cut F into Ft with two auxiliary functions ls2 and ls%\ 
The degrees of freedom concerned by Ft will have to obey simultaneously the inequalities ls2 < 0 and ls% < 0, 
while those concerned by To = F \ Ft will be the one of F which satisfy ls2 > 0 or Iss > 0. Using level-set 
functions, the position and shape of cracks can be modified when running inversions and cracks can be extended 
to study their propagation. 

To gain computation time when the crack geometry is updated several processes are followed. Let us de¬ 
note by {tpi} the standard uncut basis functions of some discrete finite element space VA C H 1 (D) introduced 
in Section 3.1. This space is made of basis functions lying in the whole domain D, which are uncut by the 
boundary F. We denote by 


( A o)ij = f^(TL(<Pi) e(<Pj)dn 

the symmetric stiffness matrix in D, which is independent of the crack and its extension. The integration method 
for computing each of these integrals is standard too. This matrix can be stored as the same time as some of 
its decompositions which enables us to invert it quickly and efficiently. From that, the basis functions {<pf} 
of the spaces (see Section 3.1 for more details) can be deduced with the help of reduction and extension 
matrices that we denote by R± and E± respectively: The reduction matrices R± enables to select the indexes 
of family {<, p i } used to define family The matrices E ± have the inverse role. Note that these matrices - 

R ± and E ± - are sparse and binary, and thus can be easily manipulated. Moreover, they satisfy the following 
useful properties: 


R+E+ m l, R + = E +T , R~E~ = I, R m E~ T . 

At this stage, we do not have to get the functions <~pf multiplied by Heaviside functions (mentioned in Section 3.1) 
as they are only used in the integrations over F. Now we define 

Aq = R^Ao, Aq = R Aq. 

Recovery of stiffness matrices A ^ in the domains (see Section 4.1) can be done from a local re-assembly 
of the integration terms for the matrices A^, by including the Heaviside functions for the terms whose basis 
functions are intersected by the level-set representing F. In practice we have easily access to the indexes which 
correspond to the local re-assemblies. 
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5 Numerical tests of convergence 

5.1 Numerical experiments without stabilization in 2D 

Rates of convergence for displacement in 2D tests 

Given a square Fl = [0; 1] x [0; 1], we consider for F a straight inclined line splitting Fl into two parts. We choose 
for T the set represented by the level-set function 

lsi(x,y) = y-2(x-x 0 ). 

The crack Ft is then chosen to be delimited by the secondary level-set functions: 


ls 2 {x,y) = x A - x, ls 2 (x,y) = x - x B , 

as the points of To satisfying ls 2 < 0 and Iss < 0. For instance we can take xq = 0.317, xa = 0.47 and 
xb = 0.52. Let us recall that we have defined To = T \ Ft- For the displacement, we consider the following 
exact solution 


U ex (x,y) 
u ex{x,y) 


( (x + y) cos(x) 
\ {x-y)sm(y) 


if lsi(x,y) > 0, 


( (x + y)cos(x) -D\{x,y) 
V (x-y)sin(y) - D 2 {x,y) 


if lsi(x, y) < 0. 


The jump D = (Di 1 D 2 ) t across Fq is no longer equal to zero. However, it is chosen to be constant, in order to 
avoid introducing jumps in normal derivatives across this boundary (see Remark 1). 

In the figures below we show the results of computation of the relative errors on the displacement, for different 
choices of the finite element spaces and W^. For instance, the couple P2/P0 indicates that we have 
considered a standard continuous P2 element for both partial displacements and u^, and a discontinuous 
P0 element for the multiplier A. We deduce an approximation of the order of convergence for the global 
displacement reconstructed afterwards from the partial ones and u^. 
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Figure 3: L 2 (fl) and H 1 (11)-relative errors (in %) on the displacement in function of the mesh size h, without 
stabilization, and estimation of convergence rates from the slope of the curves by linear regression. 


We observe (in Figure 3) that the optimal rate of convergence for the displacement seems to be better reached 
- up to round numbers errors - when PI or Q1 elements are chosen for the multiplier instead of discontinuous 
PO or QO elements, for the L 2 (Q)-relative errors as well as for the H 1 (fl)-relative errors. Furthermore, we notice 
that the rates of convergence of the H 1 (^)-norm for the displacement are better than expected when PI or 
Q1 elements are chosen for the multiplier. Last, note that some computations are missing - namely tests with 
Q3/Q1 elements and some tests with P3/P1 elements - because in these cases computing errors implies handling 
quantities which are close to the engine precision, making these tests irrelevant. 

Rates of convergence for the multiplier 

We denote by U+,, XJ~ x and A ex the discrete vectors interpolating the functions u+. := u+.| Q+ , u~ x := 
and X ex = ^<7z / (u^ c )n ± respectively. The L 2 -type error introduced by the approximation A of the exact vector 
A ex is considered as the square root of 

f |a L (U+)n+ + A| 2 dr 0 + [ |a L (U-)n- -A| 2 dr 0 . 

J r 0 J r 0 
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For practical purposes, we compute this error by developing the underlying inner product as 


(«, U + ) + (A- U U-, U- ) + 2(A+ X u+, A ex ) - 2 (A- X U- X , A ex ) + 2 <A aa A, A) (27) 

where, the matrices A± u , and A\\ are the discretized representations of the respective following bilinear 
forms 


At u : (u, v) / a L ( ujn* • cr L (v)n ± dr 0 , A± x : (u, v) / cr L (u)n ± • Adr 0 , 
Jr 0 J r 0 


-4aa 


: aO ^ / 

J r n 


A • /xdr 0 . 


Relative errors provided in Figure 4 are thus computed as the square root of the inner product (27) divided by 

11 l 2 (Fq ) 11 l 2 (Fq ) i ^ ex) "b (^uu ^ ex')^ ex)' 


Without stabilization 



Without stabilization 



Figure 4: L 2 (fl) -relative error (in %) on the multiplier in function of the mesh size h, without stabilization, and 
estimation of convergence rates from the slope of the curves by linear regression. 

In Figure 4 the optimal rate of convergence for the multiplier is obtained for PO or QO elements as well as 
for PI or Q1 elements. Note that the theoretical analysis for this fictitious domain approach does not guarantee 
a priori such a good convergence for dual variables. Anyway, without performing any stabilization, it seems 
that, for most geometric configurations the method provides us a good approximation for the values of this 
multiplier. Pathological configurations which highlight the necessity of a special treatment are in general hard 
to find. However, we show in Section 6 that the stabilization technique brings us a gain of robustness with 
respect to the geometry, as it is shown in Figure 10. 

Besides, as expected in view of Remark 1, we observed numerically that the quantity 

lki(U+)n+ + a t (U-)n-£ 2(ro) 

is negligible compared with ||<tl(U^)ii + + ||cri / (U“ a ,)n _ 11, without imposing a priori the condition 

[ctl(u)] n + across Tq. 


Illustration in 2D: Deformation in volcanic rift zones 

In order to illustrate the method on realistic 2D tests, we consider the models treated in [PDD+83], namely 
the computation of characteristics of displacements due to inside cracks. We consider the rectangular domain 
[0; 100] x [0; 50] endowed with a Cartesian mesh. The only right-hand-side is a constant pressure applied on both 
sides of the fracture as a Neumann-condition. The fracture position is determined by the points of coordinates 
(x,y) satisfying 

y = 2(x- 48.0) -1- 35.0, 35 < y < 45. 
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Standard P2 elements are chosen for the displacement, and PO elements are chosen for the multiplier. Stabiliza¬ 
tion is not performed, since here we are not interested in computing the multiplier. Illustrations are presented 
in Figure 5 and Figure 6. 



Figure 5: Representation of the intensity of displacement due to an inclined straight fracture, on a Cartesian 
mesh warped with respect to the deformation, for different numbers of subdivisions, respectively: 25 x 12, 
50 x 25, 100 x 50 and 200 x 100. 



Figure 6: Outline of level curves for the displacement due to an inclined straight fracture. 
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5.2 Choice of the stabilization parameter 7 

We are now interested in the stabilized problem (14) whose the explicit matrix formulation is given by (19). Like 
in [HR09] and [CFL14], recall that we choose 7 = 70 * h, with 70 > 0 constant, independent of the mesh size h. 
This constant has to be chosen judiciously. On one hand it represents the importance given to the quality of 
the approximation for the multiplier, but on the other hand, the stabilization term degrades the coerciveness 
of the whole system, so 70 has to remain moderate. In Figure 7 we provide a test in which the relative error on 
the multiplier A is computed for a range of values for 70 . 




Figure 7: L 2 (Fo)-relative error (in %) on A for different position of the crack Ft- 

Note in the graph on the left of Figure 7 that the first singular approximations on the multiplier occur for 
70 < 0.1, and become more frequent and chaotic for larger values of 70 . In the graph on the right we examine 
the error on A for more precise values of 70 < 0.1; We notice that the approximation of A gets better when 70 
increases, until some values of 70 > 0.04 which generate the first picks. Thus we choose 70 < 0.03 for the rest 
of the study. 

5.3 Numerical orders of convergence with stabilization 

In this subsection, we provide the same tests as performed without stabilization in Section 5.1. The stabilization 
technique is performed with the parameter 70 = 0.03. The results are given in Figure 8 and Figure 9. As expected 
for a given geometry, we observe nearly the same rates of convergence as obtained without any stabilization, 
particularly for the multiplier. 
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With stabilization 




Figure 8: L 2 (f2) and H 1 (£2)-relative errors (in %) on the displacement in function of the mesh size h, with 
stabilization, and estimation of convergence rates from the slope of the curves by linear regression. 

6 Robustness with respect to the geometry 

In order to highlight the main interest of the stabilization, we test the convergence behavior for different 
geometries, namely when the crack intersect the global mesh in different manners. For that, we perform two 
tests for which we analyze the relative error for the multiplier of L 2 (Fo) as well as relative errors on global 
displacement, in L 2 (fl) and H 1 (0), with and without using the stabilization technique. 

In the first tests we make the length of the crack vary, keeping its position, and in the second tests the length 
is kept constant, but the position is varied. In these first tests, we do not observe any significant difference 
whether we perform the stabilization or not: Results are good and similar, quite better for the multiplier with 
stabilization. These results are not shown. The second tests show the L 2 (To)-relative error for the multiplier 
with and without stabilization; They are represented in Figure 10. In these second tests the computation of 
errors on the displacement do not reveal any significant improvement due to the stabilization technique, and so 
they are not presented either. 
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With stabilization 


With stabilization 




h 


Figure 9: L 2 (fl)-relative error (in %) on the multiplier in function of the mesh size h, with stabilization, and 
estimation of convergence rates from the slope of the curves by linear regression. 




Figure 10: L 2 (Fo)-relative error (in %) on the multiplier without stabilization (in blue) and with stabilization 
(in red), for various positions of the crack, with 70 = 0.0005 (left) and 70 = 0.03 (right). 


In Figure 10, the crack determined by the points of coordinates (x, y) satisfying 

y — 2{x — xo) = 0 , xa — x < 0 , x — xb < 0 

has been considered in the square [0; 1] x [0; 1], corresponding to the exact solutions given in Section 5.1. The 
abscissas xq, xa and xb have been changed simultaneously, with xa going from 0.000 to 0.950, while keeping 
xb = xa + 0.050 and x$ = xa — 0.153. Finite elements P 2 and P0 have been chosen for the displacement and 
the multiplier respectively. 

We observe that bad computations - corresponding to huge L 2 (r 0 )-relative error - occur in fewer cases with the 
stabilization technique, and have a smaller L 2 (Fo)-relative error. Moreover, the plots show that it is delicate 
to choose the stabilization parameter 70 ; If it is chosen too large (like for instance 70 = 0.03), the problematic 
situations occur more frequently, and thus it is preferable to choose a quite small 70 (like 70 = 0.0005), even if 
the system is close to a system without any stabilization terms. 
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7 Realistic numerical simulations: Deformation of a volcano 


The code developed in 2D has been extended to dimension 3, and validated with the same type of convergence 
curves as the ones given in Section 5 by considering exact solutions. As an illustration, we consider the geo¬ 
physical problem of a pressurized crack inside a volcano. 

A 3D magma-filled crack (also called a dike) is approximated with a set of triangles, as represented in 
Figure 11. 



Figure 11: 3D magma-filled crack, called dike , represented by a set of triangles. 

Note that we do not consider this crack as a mesh: This way of representing the crack is just a practical 
means for describing its geometry from degrees of freedom - the vertices of the triangles - which are independent 
of the global mesh. 


7.1 Extension of the crack 


The purpose of this subsection is to explain how to construct an extension of the crack which splits the 3D 
computational domain, when the crack is represented by a set of triangles, such as the one given in Figure 11. 
The process of construction is the following: Given a triangle ABC whose the centroid is denoted by G, we 
define a point S whose orthogonal projection on the triangle is G. The distance between S and G is empirically 
chosen to be of the same order as the size of the triangle (namely the square root of the triangle area). We next 
define the cone of apex S based on the triangle ABC , and we cut this cone with ABC. Thus it enables us to 
define the extension of the triangle as the remaining boundary of the initial cone, and also a region delimited 
by the crack and its extension. See Figure 12. 




Figure 12: Construction of the extension of a single triangle modeling a piece of crack. 

The same process is followed for all the triangles, so that the region we define (blue color in Figure 12) is 
obtained as the reunion of all regions corresponding to each triangle. Note that the side of the triangle on which 
we have set the point S can be chosen arbitrarily for the first triangle, but for practical reason this side - in 
other words the orientation of the normal vector - has to be kept for the rest of the triangles. 
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7.2 Results and comparison with experimental data 

The mesh we use for the whole computational domain is constructed generically from a digital elevation model 
provided by IGN (Institut Geographique National, French National Geographic Institute). From these 
surface data, we first construct a surface mesh using a matlab code (Figure 13, left), and then the volume mesh 
is generated automatically with Gmsh software (Figure 13, right). 



Figure 13: The 3D mesh of the computational domain is generated from digital elevation model, from which a 
surface mesh is first generated (left) before creating the volume mesh (right). 

Note that this mesh does not take the crack into account. The local refinement we observe is due to the fact 
that some eruptive fissures corresponding to lava emission are observed at the surface (see the upper triangles 
on Figure 11), which indicates where the magma-filled crack is intersecting the ground surface. 

Elastic parameters are chosen as E = 5000 MPa for the Young modulus and v — 0.25 for the Poisson ratio 
following [CDF05]. Recall that the Lame coefficients can be related to the elastic coefficients v and E by the 
formulas below: 


Ev E 

(1 — 2i/)(1 + i/) ’ ^ 2(1 + v) 

Pressure forces are applied on both sides of the crack Ft, with values p = 5 MPa. The crater visible on the 
figures below has a diameter approximately equal to 1 km, and the crack surface is about 2.8 km 2 . The global 
mesh consists of 8376 points and 4969 tetrahedral cells. The numerical simulation which solves the displacement 
corresponding to the crack of Figure 11 in the whole computational domain is presented in Figure 14. 
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Figure 14: Amplitude of the displacement due to pressures applied on the crack. The domain is warped by 1, 
1000, 2000 and 3000 respectively, w.r.t. the displacement. The maximum - about 0.65 m - is reached at depth. 


The large deformations observed at the surface (when the domain is warped w.r.t. the displacement) are 
due to surface openings at the top of the crack (Figure 11). We can remove these surface openings, by removing 
the top triangles of the fracture of Figure 11; The corresponding displacement is represented in Figure 15. 



Figure 15: Amplitude of the displacement due to pressures applied on the fracture without surface openings. 
The domain is let unwarped (left) and then warped by 2000 (right), w.r.t. the displacement. The maximum - 
about 0.62 m - is reached in depth. 
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Our simulation can be compared qualitatively with measurements obtained by synthetic aperture radar 
[CDF05] (see Figure 16 or Figure 12 for instance), which validates qualitatively our simulation, and thus 
our method. Moreover, the order of magnitude of the displacement computed (about 0.50 m) seems also to 
correspond to the experimental observations. 

8 Conclusion 

In this work, we have proposed a new finite element method for considering a crack and simulating the dis¬ 
placement and stress induced by a crack located inside a heterogeneous and anisotropic elastic medium. In 
order to avoid re-meshing and save computation time when inverting from surface deformation data or when 
simulating a propagating fracture, we consider a single global mesh and we develop a fictitious domain approach. 
Numerical tests were performed in order to highlight the optimal accuracy of the method. Robustness with 
respect to the geometry is ensured by the implementation of a stabilization technique. A theoretical analysis 
underlines the relevance of the latter for computing a good approximation of a dual variable. Numerical tests 
enable us to think that our method is practical for an algorithmic framework in which multiple ways of taking 
into account the crack will be considered. Indeed, the inversion of the position/shape of the crack - inside a 
volcano, from surface measurements for instance - is an algorithmic problem that can involve this kind of dual 
variables (shape derivatives techniques involving adjoint direct systems for instance). The illustration of this 
method on an inverse problem will be presented in a forthcoming work. The fictitious domain approach we have 
developed enables us to consider a single global mesh for all kind of geometries associated with the crack. In 
the context of inverse problems, we expect its simplicity to be quite effective in terms of computational time. 
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